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ABSTRACT 

The goal of this research is to investigate how well various turbulence models 
can describe physical properties of the upper convective boundary layer of the 
Sun. An accurate modeling of the turbulence motions is necessary for understand- 
ing the excitation mechanisms of solar oscillation modes. We have carried out 
realistic numerical simulations using several different physical Large Eddy Simu- 
lation (LES) models (Hyperviscosity approach, Smagorinsky, and dynamic mod- 
els) to investigate how the differences in turbulence modeling affect the damping 
and excitation of the oscillations and their spectral properties and compare with 
observations. We have first calculated the oscillation power spectra of radial and 
non-radial modes supported by the computational box with the different turbu- 
lence models. Then we have calculated the work integral input to the modes 
to estimate the influence of the turbulence model on the depth and strength of 
the oscillation sources. We have compared these results with previous studies 
and with the observed properties of solar oscillations. We find that the dynamic 
turbulence model provides the best agreement with the helioseismic observations. 

Subject headings: convection — methods: numerical — Sun: oscillations 
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1. Introduction 

Dominant acoustic sources within the Sun are generated by strong fluctuations in the 
outer convective layers. Turbulent motions stochastically excite the resonant modes via 
Reynolds stresses and entropy fluctuations. The dominant driving comes from the interac- 
tion o f the nonadiabatic, incoh erent pressure fluctuations with the coherent mode displace- 



ment (INordlund fc Stein! 120011 ) . The modes excitation sources occur cl ose to the surface. 



main ly in the intergranular lanes and near the boundaries of granules (IStein fc Nordlund 



20011 ) . Thus an accurate modeling of the turbulence motions is necessary to understand the 
excitation mechanisms of solar oscillation modes. The correct choice of turbulence model is 
also important in many other astrophysical simulations. 

The objective of this research is to study the influence of turbulence models on the exci- 
tation mechanisms by means of realistic numerical simulations. We have compared different 
physical Large Eddy Simulation (LES) models (Hyperviscosity approach, Smagorinsky, and 
dynamic models) to show the influence on the damping and excitation of the oscillations. 

The organization of this paper is as follows. In §2, we describe the main lines of the 
code and the different turbulence models. The kinetic energy of the radial modes obtained 
with the different turbulence models are presented in §3. Then a comparison of the results 
obtained with the different turbulence models for the non-radial modes is given in §4. The 
work integral input to the modes is calculated in §5 in order to estimate the influence of the 
turbulence models on the depth of the oscillation sources. 



2. Numerical model 

We use a 3D, compressible, non-linear radiative-hydrodynamics code developed by Dr. 
A. Wray for simulating the upper solar photosphere and lower atmosphere. This code 
takes into account several physical phenomena: compressible fluid flow in a highly strat- 
ified medium, radiative energy transfer between the fluid elements, and a real-gas equation 
of state. The equations we solve are the grid-cell average (henceforth called "average") 
conservation of mass, 

dp 



conservation of momentum, 

dpuj 
dt 



, )t + (p Ui ). = 0, (1) 



+ (pUiUj + (Pij + pTij)) d = -p<f>,i, (2) 
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and conservation of energy, 
dE 

— + {Em + (Pij + pr^uj - (k + K T )Ti + F[ ad ) . = 0, (3) 

where p denotes the average mass density, Ui the Favre-averaged velocity, and E the average 
total energy density E = \pu{iLi + pe + p<p, where is the gravitational potential and e 
is the Favre-averaged internal energy density per unit mass. F™ d is the radiative flux, 
calculated by solution of the radiative transfer equation, and is the average stress tensor 
Pij = (p + 2/xWfc,A:/3) Sij — p (ui t j + Ujj) with p the viscosity. The fluid pressure p is a function 
of e and p through a tabulated equation of state; is the Reynolds stress, k is the molecular 
thermal conductivity, and kt is the turbulent thermal conductivity. 



The turbulence mod els used are the origin al Smagorinsky model (ISmagorinskylll963l ) and 



its dynamic formulation (IGermano et al. Il99lh . herein called simply the dynamic model. The 



Reynolds stresses Tjj are modeled as in the usual Samgorinsky formulation by writing them 
in terms of an eddy viscosity formed from the large-scale stress tensor Sij = (uij + Uj^)/2: 

n, = -2C s tf\S\(S VJ - u k , k Sij/3) + 2C c A 2 \S\%/3, (4) 

where \S\ = ^2SijSij and A = (AxAyAz) 1 ^ 3 , Ax, Ay, and Az are the grid step sizes. 
The two parameters, C$ and Cc, which are the classical Smagorinsky coefficient as used in 
incompressible flow and a coefficient associated with the trace of the subgrid Reynolds stress 
(which is absent in incompressible flow), must be specified in some way. Constant values 
were used in some runs and in others these parameters were determined by the dynamic 
method using planar averages. The turbulent Prandtl number was taken as unity to set Kt- 
The molecular viscosity p and thermal conductivity k were taken to be zero as their solar 
values are exceedingly small. 

We simulate the upper layers of the convection zone using 66x66x40 grid cells. The 
region extends 6x6 Mm horizontally and from 2.5 Mm below the visible surface to 0.5 Mm 
above the surface. This co mputational box have bee n chosen to directly compare with the 



previous results obtained by IStein fc Nordfundl (120011 ) who used a numerical viscosity model 



not related to a particular turbulence model. 



3. Kinetic energy of radial modes 

First we have studied how the kinetic energy is dissipated for the different turbulence 
models. Specifically we have calculated the oscillation power spectra of radial modes. Those 
modes are extracted by horizontal averaging of the vertical velocity and Fourier transform in 
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time. The results presented here (Figured] and Figure[2]) have been obtained with simulations 
of 60 hours of solar time, and the snapshots were saved every 30 seconds. 

Three modes can be clearly seen in the spectra of the horizontally averaged, depth- 
integrated kinetic energy obtained with all the three turbulence models (Figure [2] (left 
panel)), they correspond to the maxima in the kinetic energy. The resonant frequencies 
supported by the computa tional box are 2.6, 4.0, 5 .6 mHz. These frequencies are very close 



to the values obtained by lStein fc Nordlundl (120011 ). We see that for the hyperviscosity ap 



proach and the dynamic model the modes within the computation box are excited with the 
same magnitude whereas the Smagorinsky model yields a lower magnitude. 

The kinetic energy spectra as a function of frequency and depth (Figure [1]) confirm that 
the dissipation is weaker with the minimal hyperviscosity approach. In this case, almost 
nearly the numerical dissipation plays a role. Thus the kinetic energy is higher for high 
frequencies in comparison with the results obtained with the other turbulence models. The 
spectrum obtained with the enhanced hyperviscosity approach and the dynamic model show 
a weakly higher dissipation compared to the calculation with the minimal hyperviscosity 
approach. Moreover the dissipation applied by the hyperviscosity approach and the dynamic 
model does not interact with the three oscillations modes. In fact, the dissipation scale is 
smaller than the scale of the acoustic modes. With the Smagorinsky model, the excitation 
of the modes is weaker and the scale of the dissipation is close to the scale of the third mode. 
The patterns do not extend above 6 mHz (against 12 mHz with the hyperviscosity approach 
and 10 mHz with the dynamic model). The cut-off frequency obta ined with the Smagorinsky 



model is very similar of that obtained by lStein fc Nordlundl (120011 ): the kinetic energy is very 



low above 6 mHz. Their calculations have been performed with artificial numerical viscosity. 

These results show that the spectral kinetic energy of the radial modes is rather sim- 
ilar for the hyperviscosity and dynamic model, but the Smagorinsky model is too much 
dissipative. 



4. Kinetic energy of non-radial modes 

The non-radial modes have been extracted by performing 2d spatial Fourier transform 
on the surfaces of vertical velocity at each time step. Thus we obtain a power spectrum 
for a particular horizontal wave number k\ = k\ + k 2 y . Then we take Fourier transform in 
time. We have especially considered modes with horizontal wavelength corresponding to 
the box size L = 6 Mm. That corresponds to angular degree of I ~ 740 ikh = 1 MvrC x \ 
Figure [2] (right panel) shows the power spectra for the different turbulent models. The 
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leftmost peak corresponds to the surface gravity (f) mode, while the others correspond to 
acoustic (p) modes. The f-mode peak has the same frequency independent of the depth of 
the computational box (2.8 mHz). One can see that the Smagorinsky model gives lowest 
power, and the modes are less excited. The hyperviscosity approach and the dynamic model 
give very close results. The modes have approximately the same magnitude. The influence 
of the turbulence models is similar for the radial and non-radial modes. 



5. Calculation of the p-mode excitation rates 



The mode excitation rate is calculated using the same method presented by lStein &: Nordlund 



(120011 ) . The rate of energy input to the modes per unit surface area (erg. cm 2 .s 1 ) is 

A < E u > _ u 2 \ j r dr5P: jdUdr) \ 2 
At 8AzaEL 



(5) 



where 5P* is the Fourier transform of the nonadiabatic total pressure. Av the frequency in- 
terval for the Fourier transform. £ w is the mode displacement f or the radial mode of angular 



freque ncy uj. It is obtained from the eigenmode calculations of IChristensen-Dalsgaard et al 



( 119961 ). His spherically symmetric model S gives 35 radial modes what provide much bet- 
ter frequency coverage in comparison with the three resonant modes obtained within the 
simulation box. E u is the mode energy per unit surface area (erg.cm~ 2 ) defined as: 

1 .2 /\ , 2(r\ 2 



E U = -(S J drpZJ{-) . (6) 

We have calculated the rate of stochastic energy input to modes for the entire solar 
surface in order to compare with the observed results (Figure E]). The rate of energy input 
to the solar modes is obtained by multiplying the rate of the simulation modes by the ratio 
of the mode mass of the solar modes over the mode mass in the computational domain. 
The comparison between observed and calculated rate shows that the Smagorinsky model 
gives too low values (more than one decade below the observations). When we consider the 
minimal hyperviscosity approach the main discrepancy is located between 2 and 3 mHz. We 
can see a peak in this range, which is much less pronounced in the observed results. This 
peak is present also with the enhanced hyperviscosity approach but it is lower. It is no 
longer present with the dynamic model. The best agreement is obtained below 3 mHz using 
the dynamic model. Above 3 mHz, both the enhanced hyperviscosity approach and dynamic 
model provide good results. The origin of the peak o r plateau in th e sol ar power spectrum in 



terms of the mixing length theory was discussed by iGoughl (Il980l ) and iBalmforth &: Gough 



(Il990h 
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The distributions of the integrand of the work integral as a function of depth and 
frequency (Figure H|) c an explain the pres e nce o f this peak. The distributions are similar to 
the results obtained by lStein fc Nordlundl (120011 ). The most driving is concentrated between 
the surface and 500 km depth around 3-4 mHz. We can see that the main differences between 
the magnitude obtained with the different turbulence models are located in the region around 
2 Mm. We can observe that the excitation obtained with the minimal hyperviscosity (Figure 
S] (top left panel)) is high between 2 and 3 mHz at each vertical position. Conversely the 
excitation decreases as the depth increases with the enhanced hyperviscosity approach and 
the dynamic model. The excitation magnitude becomes low when the depth is greater than 
2 Mm with the enhanced hyperviscosity approach (1.5 Mm with the dynamic model). These 
differences explain the difference in magnitude of the peak between 2 and 3 mHz for the 
rate of stochastic energy. The presence of the peak in the observations underlines that the 
excitation is only located close to the visible surface of the Sun. 



6. Conclusion 



The goals of this research was to investigate how well various turbulence models can 
describe the convective properties of the upper boundary layer of the Sun and to study the 
excitation and damping of acoustic oscillations. Results obtained with the hyperviscosity 
approach have been compared with those obtained with the Smagorinsky and dynamic tur- 
bulence models. We have seen that the dissipation is very high with the Smagorinsky model 
while the hyperviscosity approach and dynamic modes give similar results. Besides we find 
that the dynamic turbulence model provides the best agreement with observations. 
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Fig. 1. — Logarithm of kinetic energy as a function of depth and frequency (erg.cm~ 3 ). 
Top left: Minimal hyperviscosity approach. Top right: Enhanced hyperviscosity approach. 
Bottom left: Smagorinsky model (Cs=0.2). Bottom right: Dynamic model. 
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Fig. 2. — Left: Logarithm of the kinetic energy, horizontally averaged and integrated over 
depth for the different turbulence models I = (erg.cm~ 2 ). Right: Power spectra density of 
vertical velocity with different turbulent modes for angular degree I = 740 at 0.5 Mm below 
the surface ((cm.s -1 ) 2 ). 
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Fig. 3. — Comparison of observed and calculated rate of stochastic energy input to modes 
for the entire solar s urface (erq.s~ l ). Observ ed distributions (squares) come from SOHO 
GOLF for 1 = - 3 faoca Cortes et adll99flh 
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Fig. 4. — Logarithm of the work integrand (eq. (jSJ) in units of erg.cm~ 2 .s" 1 ), as a function 
of depth and frequency. Top left: Minimal hyperviscosity approach. Top right: Enhanced 
hyperviscosity approach (A=0. 4). Bottom left: Smagorinsky model ((7,5=0.2). Bottom right: 
Dynamic model. 



